options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

hierarchical model

class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)
n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
  scale_shape_manual(values=1:k)

qplot(x,y,col=c)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

estimate as no class

ex8-0.stan

y~N(b00+b10*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-0.stan') 
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)
## Initial log joint probability = -199299 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       38      -156.345     0.0262739    0.00171095           1           1       86    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__    -156.34
##    b0         6.78
##    b1         1.84
##    s         13.83
##    m0[1]     35.77
##    m0[2]     29.30
##    m0[3]     21.80
##    m0[4]      8.04
##    m0[5]     43.23
##    m0[6]     37.14
##    m0[7]      8.68
##    m0[8]     32.52
##    m0[9]     36.62
##    m0[10]    36.72
##    m0[11]    18.70
##    m0[12]    42.85
##    m0[13]    43.41
##    m0[14]    25.13
##    m0[15]    12.52
##    m0[16]    18.77
##    m0[17]    11.51
##    m0[18]    31.43
##    m0[19]     7.65
##    m0[20]    32.54
##    m0[21]    31.77
##    m0[22]    28.54
##    m0[23]    42.99
##    m0[24]    29.98
##    m0[25]    43.39
##    m0[26]    33.41
##    m0[27]    42.39
##    m0[28]    41.15
##    m0[29]    13.95
##    m0[30]     8.30
##    m0[31]    25.88
##    m0[32]    39.08
##    m0[33]    40.99
##    m0[34]    14.54
##    m0[35]    21.13
##    m0[36]     9.55
##    m0[37]    30.21
##    m0[38]    23.15
##    m0[39]    26.80
##    m0[40]     7.31
##    m0[41]     9.39
##    m0[42]    38.97
##    m0[43]    19.64
##    m0[44]    22.56
##    m0[45]    35.67
##    m0[46]    34.55
##    m0[47]    34.54
##    m0[48]    34.14
##    m0[49]    25.23
##    m0[50]    35.19
##    y0[1]     35.35
##    y0[2]     23.41
##    y0[3]     23.58
##    y0[4]     -2.21
##    y0[5]     45.05
##    y0[6]     29.95
##    y0[7]    -17.55
##    y0[8]     17.57
##    y0[9]     26.08
##    y0[10]    33.38
##    y0[11]     7.80
##    y0[12]    36.52
##    y0[13]    29.11
##    y0[14]     5.03
##    y0[15]     0.58
##    y0[16]    -8.45
##    y0[17]    13.32
##    y0[18]    23.72
##    y0[19]    17.24
##    y0[20]    41.35
##    y0[21]     4.75
##    y0[22]    21.90
##    y0[23]    46.17
##    y0[24]    48.29
##    y0[25]    16.93
##    y0[26]    13.42
##    y0[27]    45.88
##    y0[28]    31.25
##    y0[29]    15.02
##    y0[30]    15.62
##    y0[31]    19.26
##    y0[32]    35.38
##    y0[33]    37.60
##    y0[34]    24.29
##    y0[35]    26.81
##    y0[36]   -15.65
##    y0[37]    44.23
##    y0[38]    25.64
##    y0[39]     9.98
##    y0[40]   -13.88
##    y0[41]    -3.61
##    y0[42]    45.82
##    y0[43]     5.23
##    y0[44]    19.56
##    y0[45]    51.28
##    y0[46]    55.60
##    y0[47]    18.45
##    y0[48]    38.34
##    y0[49]    11.28
##    y0[50]    44.00
##    m1[1]      7.31
##    m1[2]     11.32
##    m1[3]     15.33
##    m1[4]     19.34
##    m1[5]     23.35
##    m1[6]     27.36
##    m1[7]     31.37
##    m1[8]     35.39
##    m1[9]     39.40
##    m1[10]    43.41
##    y1[1]      6.80
##    y1[2]     22.67
##    y1[3]     21.33
##    y1[4]     47.35
##    y1[5]     30.68
##    y1[6]     15.24
##    y1[7]     14.36
##    y1[8]     21.25
##    y1[9]     33.90
##    y1[10]    65.01
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -155.23 -154.91  1.27  0.96 -157.63 -153.89 1.01      487      470
##    b0        6.72    6.74  4.23  4.05    0.24   13.63 1.03      221      275
##    b1        1.84    1.84  0.33  0.32    1.30    2.38 1.02      268      274
##    s        14.45   14.34  1.51  1.49   12.14   17.12 1.00     2133     1398
##    m0[1]    35.72   35.72  2.51  2.54   31.61   39.69 1.00     1542     1184
##    m0[2]    29.25   29.21  2.05  2.03   25.94   32.52 1.00     1149      980
##    m0[3]    21.75   21.76  2.26  2.24   18.17   25.26 1.01      314      364
##    m0[4]     7.98    8.00  4.04  3.88    1.74   14.59 1.03      222      274
##    m0[5]    43.18   43.16  3.48  3.43   37.45   48.86 1.00      717      955
##    m0[6]    37.09   37.12  2.67  2.65   32.71   41.30 1.00     1439     1166
##    m0[7]     8.62    8.66  3.94  3.77    2.48   15.06 1.03      222      277
##    m0[8]    32.47   32.46  2.22  2.22   28.87   36.06 1.00     1599     1337
##    m0[9]    36.57   36.60  2.61  2.61   32.30   40.70 1.00     1495     1181
##    m0[10]   36.67   36.70  2.62  2.62   32.38   40.82 1.00     1486     1181
##    m0[11]   18.65   18.66  2.56  2.54   14.58   22.73 1.02      261      315
##    m0[12]   42.80   42.77  3.43  3.38   37.15   48.36 1.00      739      996
##    m0[13]   43.36   43.32  3.51  3.45   37.58   49.07 1.00      708      936
##    m0[14]   25.08   25.13  2.06  2.12   21.79   28.36 1.00      457      448
##    m0[15]   12.47   12.51  3.36  3.23    7.17   17.83 1.03      227      290
##    m0[16]   18.71   18.73  2.55  2.53   14.66   22.78 1.02      261      315
##    m0[17]   11.46   11.52  3.51  3.35    5.97   17.09 1.03      225      290
##    m0[18]   31.38   31.37  2.14  2.16   27.91   34.85 1.00     1560     1364
##    m0[19]    7.60    7.62  4.10  3.92    1.30   14.31 1.03      221      275
##    m0[20]   32.49   32.47  2.22  2.21   28.88   36.08 1.00     1600     1337
##    m0[21]   31.71   31.70  2.16  2.18   28.24   35.20 1.00     1575     1367
##    m0[22]   28.49   28.44  2.03  2.03   25.20   31.73 1.00      896      951
##    m0[23]   42.94   42.91  3.45  3.40   37.26   48.56 1.00      730      955
##    m0[24]   29.92   29.92  2.07  2.06   26.56   33.25 1.00     1365     1275
##    m0[25]   43.34   43.31  3.51  3.45   37.56   49.05 1.00      709      936
##    m0[26]   33.35   33.35  2.29  2.28   29.62   37.07 1.00     1606     1280
##    m0[27]   42.34   42.30  3.36  3.32   36.79   47.79 1.00      769      965
##    m0[28]   41.10   41.06  3.19  3.14   35.81   46.25 1.00      865     1018
##    m0[29]   13.90   13.93  3.16  3.06    8.90   18.92 1.03      231      292
##    m0[30]    8.24    8.29  4.00  3.86    2.05   14.78 1.03      222      274
##    m0[31]   25.83   25.89  2.04  2.11   22.58   29.08 1.00      514      518
##    m0[32]   39.03   39.01  2.91  2.91   34.17   43.62 1.00     1135     1040
##    m0[33]   40.95   40.89  3.16  3.12   35.68   46.05 1.00      881     1018
##    m0[34]   14.48   14.50  3.08  2.98    9.60   19.39 1.03      233      299
##    m0[35]   21.07   21.08  2.32  2.29   17.36   24.71 1.01      299      342
##    m0[36]    9.49    9.52  3.80  3.62    3.54   15.73 1.03      223      283
##    m0[37]   30.16   30.16  2.08  2.08   26.80   33.53 1.00     1415     1258
##    m0[38]   23.10   23.15  2.16  2.17   19.68   26.48 1.00      356      388
##    m0[39]   26.74   26.77  2.02  2.07   23.51   29.99 1.00      605      651
##    m0[40]    7.25    7.27  4.15  3.96    0.89   14.02 1.03      221      275
##    m0[41]    9.34    9.35  3.83  3.63    3.34   15.61 1.03      223      283
##    m0[42]   38.92   38.92  2.89  2.89   34.09   43.49 1.00     1153     1040
##    m0[43]   19.58   19.60  2.46  2.44   15.65   23.46 1.02      273      334
##    m0[44]   22.51   22.54  2.20  2.16   19.01   25.93 1.01      336      378
##    m0[45]   35.62   35.62  2.50  2.52   31.53   39.58 1.00     1546     1184
##    m0[46]   34.50   34.51  2.39  2.37   30.57   38.32 1.00     1585     1296
##    m0[47]   34.49   34.50  2.39  2.37   30.56   38.31 1.00     1585     1296
##    m0[48]   34.09   34.09  2.35  2.33   30.24   37.88 1.00     1595     1278
##    m0[49]   25.18   25.23  2.06  2.11   21.90   28.46 1.00      464      449
##    m0[50]   35.14   35.15  2.45  2.45   31.14   39.01 1.00     1563     1200
##    y0[1]    34.96   35.27 14.93 14.19    9.77   59.35 1.00     2063     1972
##    y0[2]    29.36   29.27 14.74 14.45    5.57   53.34 1.00     2055     1769
##    y0[3]    21.39   21.26 14.31 14.01   -1.60   45.12 1.00     1840     1905
##    y0[4]     7.72    7.74 14.69 14.78  -16.94   31.71 1.00     1789     1731
##    y0[5]    43.17   43.42 14.59 14.04   18.97   66.35 1.00     1694     1704
##    y0[6]    37.16   37.39 15.03 14.26   12.08   62.21 1.00     1824     1896
##    y0[7]     8.08    8.03 15.07 15.30  -16.37   32.22 1.00     1213     1683
##    y0[8]    32.22   32.37 14.31 13.92    8.73   56.76 1.00     2053     1973
##    y0[9]    36.59   37.03 14.69 14.09   12.26   61.00 1.00     1971     1818
##    y0[10]   36.44   36.30 14.57 14.42   12.49   60.49 1.00     2011     2012
##    y0[11]   18.92   19.04 14.65 14.13   -5.06   43.58 1.00     1675     1844
##    y0[12]   42.43   42.38 15.45 15.16   16.47   67.25 1.00     1844     1847
##    y0[13]   43.29   42.90 14.66 14.90   20.46   68.24 1.00     1932     2030
##    y0[14]   25.59   25.59 14.69 14.24    1.49   49.29 1.00     1843     1905
##    y0[15]   11.71   11.94 15.29 15.18  -13.78   36.49 1.00     1648     1879
##    y0[16]   18.60   18.91 14.84 14.60   -5.70   42.02 1.00     1925     1683
##    y0[17]   11.28   10.94 15.08 14.58  -13.19   36.95 1.00     1745     1839
##    y0[18]   31.31   31.33 14.32 14.34    8.31   55.36 1.00     2056     2006
##    y0[19]    8.00    7.99 14.93 14.76  -17.14   33.41 1.00     1034     1705
##    y0[20]   32.24   32.20 14.77 14.42    8.33   56.04 1.00     1751     1754
##    y0[21]   31.83   31.43 14.72 14.32    7.28   55.90 1.00     1953     1995
##    y0[22]   28.43   28.39 14.24 13.99    4.87   51.38 1.00     2016     1934
##    y0[23]   43.22   43.49 15.15 14.72   17.60   68.09 1.00     1989     1722
##    y0[24]   30.75   30.99 14.71 15.19    6.34   54.50 1.00     1900     2015
##    y0[25]   43.02   42.69 15.14 15.02   18.45   68.81 1.00     2072     1875
##    y0[26]   33.70   33.65 14.75 14.89   10.40   57.29 1.00     1729     1636
##    y0[27]   42.30   42.26 15.36 15.41   16.96   67.53 1.00     2081     1927
##    y0[28]   40.85   40.86 15.00 14.69   16.20   65.92 1.00     1869     1757
##    y0[29]   14.16   14.16 14.66 14.97   -9.57   39.06 1.00     1640     1910
##    y0[30]    8.11    8.26 14.91 14.10  -16.39   33.11 1.00     1476     1577
##    y0[31]   25.65   25.86 14.56 14.18    0.96   50.46 1.00     1768     1968
##    y0[32]   38.73   39.04 15.43 15.24   13.36   64.20 1.00     1892     2040
##    y0[33]   41.40   41.18 14.56 14.27   17.39   65.04 1.00     1958     1711
##    y0[34]   14.49   14.03 14.80 15.00   -9.14   39.36 1.00     1714     1916
##    y0[35]   21.43   21.23 14.43 14.60   -1.78   45.22 1.00     1775     1930
##    y0[36]    9.58    9.75 14.87 14.59  -15.31   34.08 1.00     1481     1891
##    y0[37]   30.15   30.34 14.56 14.47    6.20   53.72 1.00     2023     1879
##    y0[38]   23.49   23.62 14.72 14.60   -0.71   47.52 1.00     1904     2056
##    y0[39]   26.73   26.83 14.77 14.81    2.84   51.30 1.00     2122     1928
##    y0[40]    7.48    7.29 15.18 14.82  -17.19   32.84 1.00     1708     1853
##    y0[41]    9.83    9.60 14.86 15.40  -14.98   33.65 1.00     1725     1970
##    y0[42]   39.37   39.32 14.61 14.13   14.33   64.21 1.00     1798     1932
##    y0[43]   19.27   19.17 14.86 14.23   -5.61   42.94 1.00     1676     1971
##    y0[44]   22.47   22.54 14.75 14.27   -2.39   46.70 1.00     1997     1880
##    y0[45]   35.85   36.18 14.60 14.72   11.69   59.78 1.00     2004     2028
##    y0[46]   33.96   34.18 14.72 13.98    9.09   57.90 1.00     1846     2013
##    y0[47]   34.48   34.50 14.43 14.13   10.44   57.66 1.00     1957     1952
##    y0[48]   33.84   33.63 14.97 14.54    9.09   59.01 1.00     1877     1931
##    y0[49]   25.36   25.26 14.42 14.42    1.58   49.08 1.00     1792     1765
##    y0[50]   34.96   34.82 14.59 14.02   11.31   59.34 1.00     2067     2196
##    m1[1]     7.25    7.27  4.15  3.96    0.89   14.02 1.03      221      275
##    m1[2]    11.26   11.32  3.54  3.38    5.74   16.96 1.03      225      294
##    m1[3]    15.27   15.30  2.97  2.86   10.55   20.02 1.03      236      305
##    m1[4]    19.29   19.30  2.49  2.46   15.31   23.21 1.02      268      332
##    m1[5]    23.30   23.34  2.15  2.15   19.91   26.67 1.00      363      388
##    m1[6]    27.31   27.33  2.02  2.03   24.11   30.56 1.00      678      757
##    m1[7]    31.32   31.31  2.14  2.15   27.86   34.79 1.00     1557     1384
##    m1[8]    35.34   35.33  2.47  2.48   31.31   39.24 1.00     1557     1185
##    m1[9]    39.35   39.32  2.95  2.91   34.40   44.00 1.00     1082     1076
##    m1[10]   43.36   43.32  3.51  3.45   37.58   49.07 1.00      708      936
##    y1[1]     6.92    7.19 15.27 15.12  -18.30   31.01 1.00     1204     1764
##    y1[2]    10.97   10.99 15.28 14.66  -14.42   36.83 1.00     1168     1708
##    y1[3]    15.26   15.51 14.70 14.40   -8.92   38.48 1.00     1360     1595
##    y1[4]    19.31   19.33 14.61 13.71   -4.58   43.03 1.00     1684     1822
##    y1[5]    23.39   23.62 14.21 13.76   -0.28   46.91 1.00     2110     1920
##    y1[6]    27.70   27.36 14.54 13.69    3.73   51.91 1.00     2065     1783
##    y1[7]    31.84   32.06 14.92 14.88    7.97   56.48 1.00     2099     1857
##    y1[8]    35.15   35.27 14.57 14.38   10.76   58.66 1.00     2019     1785
##    y1[9]    39.38   39.16 15.30 14.82   14.23   64.87 1.00     1907     1870
##    y1[10]   43.46   43.33 14.89 14.76   18.93   67.66 1.00     2029     1988

estimate as independent class

all b0l,b1l do not have a distribution  
class c=1~k 

ex12-1.stan

yc~N(b0c+b1c*x,s)

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  vector[K] b0;
  vector[K] b1;
  real<lower=0> s;
}
model{
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-1.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -104749 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       99      -87.4726      0.128499      0.386209           1           1      151    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199      -87.4689    0.00239103     0.0123312      0.3276           1      264    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      201      -87.4689    0.00222475    0.00169226           1           1      266    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -87.47
##    b0[1]      9.92
##    b0[2]      4.67
##    b0[3]     13.60
##    b0[4]     -9.46
##    b0[5]     13.97
##    b0[6]     15.62
##    b0[7]      7.27
##    b0[8]      9.05
##    b0[9]     12.15
##    b1[1]      2.81
##    b1[2]      0.95
##    b1[3]      2.14
##    b1[4]      2.02
##    b1[5]      2.70
##    b1[6]      2.13
##    b1[7]      1.57
##    b1[8]     -0.38
##    b1[9]      0.91
##    s          3.49
##    m0[1]     32.04
##    m0[2]     39.81
##    m0[3]      5.96
##    m0[4]      8.35
##    m0[5]     57.91
##    m0[6]      2.81
##    m0[7]     16.76
##    m0[8]     24.86
##    m0[9]     26.89
##    m0[10]    48.45
##    m0[11]     3.65
##    m0[12]    30.19
##    m0[13]    65.86
##    m0[14]    36.91
##    m0[15]     7.87
##    m0[16]    17.51
##    m0[17]    14.49
##    m0[18]    24.33
##    m0[19]    11.25
##    m0[20]    43.58
##    m0[21]    24.49
##    m0[22]    22.90
##    m0[23]    30.35
##    m0[24]    16.04
##    m0[25]    58.09
##    m0[26]    44.59
##    m0[27]    56.93
##    m0[28]    22.38
##    m0[29]    21.95
##    m0[30]    15.37
##    m0[31]    35.83
##    m0[32]    61.39
##    m0[33]    22.30
##    m0[34]     8.67
##    m0[35]    19.23
##    m0[36]    16.82
##    m0[37]    27.29
##    m0[38]    13.11
##    m0[39]    22.04
##    m0[40]    -8.87
##    m0[41]    18.66
##    m0[42]    25.93
##    m0[43]    30.54
##    m0[44]    19.94
##    m0[45]    56.39
##    m0[46]    45.92
##    m0[47]    52.31
##    m0[48]    54.15
##    m0[49]    21.26
##    m0[50]     3.21
##    y0[1]     31.34
##    y0[2]     42.24
##    y0[3]      4.15
##    y0[4]     11.03
##    y0[5]     57.74
##    y0[6]      2.90
##    y0[7]     18.33
##    y0[8]     24.98
##    y0[9]     24.09
##    y0[10]    45.35
##    y0[11]     0.99
##    y0[12]    33.97
##    y0[13]    72.31
##    y0[14]    36.21
##    y0[15]    10.63
##    y0[16]    16.11
##    y0[17]    15.55
##    y0[18]    20.10
##    y0[19]    12.28
##    y0[20]    41.61
##    y0[21]    24.10
##    y0[22]    21.15
##    y0[23]    32.76
##    y0[24]    18.54
##    y0[25]    58.97
##    y0[26]    38.41
##    y0[27]    59.37
##    y0[28]    20.33
##    y0[29]    22.44
##    y0[30]    19.07
##    y0[31]    32.60
##    y0[32]    59.67
##    y0[33]    24.00
##    y0[34]    11.27
##    y0[35]    17.70
##    y0[36]    21.32
##    y0[37]    27.99
##    y0[38]     8.81
##    y0[39]    25.00
##    y0[40]   -13.42
##    y0[41]    21.28
##    y0[42]    30.65
##    y0[43]    32.41
##    y0[44]    15.24
##    y0[45]    55.36
##    y0[46]    50.97
##    y0[47]    54.41
##    y0[48]    55.49
##    y0[49]    18.54
##    y0[50]     0.42
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 1.5 seconds.
## Chain 4 finished in 1.4 seconds.
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 1.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.5 seconds.
## Total execution time: 1.8 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad      q5    q95 rhat ess_bulk ess_tail
##    lp__   -98.10 -97.65 3.89 3.74 -105.02 -92.64 1.00      623     1086
##    b0[1]    9.93   9.87 4.62 4.70    2.54  17.27 1.00     1925     1389
##    b0[2]    4.76   4.75 5.00 4.85   -3.24  12.98 1.00     1602     1272
##    b0[3]   13.63  13.68 2.99 2.94    8.87  18.54 1.00     2880     1570
##    b0[4]   -9.48  -9.39 3.80 3.68  -15.76  -3.18 1.00     2324     1476
##    b0[5]   14.13  14.11 4.90 4.98    5.92  22.38 1.00     1801     1257
##    b0[6]   15.62  15.65 3.73 3.50    9.40  21.79 1.00     2324     1364
##    b0[7]    7.30   7.31 4.12 4.16    0.66  13.93 1.00     1923     1104
##    b0[8]    9.02   9.09 5.24 5.37    0.22  17.24 1.00     1407     1275
##    b0[9]   12.19  12.08 4.29 4.36    5.44  19.30 1.00     1951     1363
##    b1[1]    2.82   2.82 0.32 0.31    2.29   3.34 1.00     1946     1221
##    b1[2]    0.94   0.95 0.35 0.35    0.35   1.51 1.00     1677     1597
##    b1[3]    2.14   2.13 0.26 0.26    1.71   2.57 1.00     2581     1653
##    b1[4]    2.02   2.02 0.26 0.24    1.60   2.44 1.00     2237     1306
##    b1[5]    2.68   2.68 0.36 0.36    2.10   3.26 1.00     1850     1398
##    b1[6]    2.13   2.12 0.26 0.24    1.72   2.55 1.00     2389     1685
##    b1[7]    1.57   1.57 0.39 0.39    0.93   2.19 1.00     1951     1477
##    b1[8]   -0.38  -0.38 0.43 0.42   -1.05   0.35 1.01     1566     1441
##    b1[9]    0.90   0.91 0.38 0.38    0.28   1.52 1.00     1949     1617
##    s        4.53   4.46 0.60 0.57    3.65   5.63 1.00     1078     1349
##    m0[1]   32.10  32.10 3.42 3.49   26.59  37.69 1.00     1910     1704
##    m0[2]   39.83  39.84 1.61 1.58   37.21  42.43 1.00     2191     1232
##    m0[3]    5.93   6.06 2.59 2.58    1.65  10.03 1.00     1389     1490
##    m0[4]    8.38   8.40 3.90 3.93    2.11  14.70 1.00     1916     1102
##    m0[5]   57.86  57.79 2.65 2.52   53.56  62.29 1.00     2124     1283
##    m0[6]    2.79   2.71 3.24 3.17   -2.36   8.17 1.00     1958     1365
##    m0[7]   16.90  16.92 4.58 4.66    9.32  24.62 1.00     1804     1258
##    m0[8]   24.84  24.86 1.92 1.79   21.75  28.00 1.00     2028     1611
##    m0[9]   26.85  26.83 2.57 2.43   22.66  30.94 1.00     1996     1585
##    m0[10]  48.47  48.49 2.26 2.24   44.76  52.09 1.00     2195     1357
##    m0[11]   3.62   3.66 2.47 2.41   -0.50   7.58 1.00     2256     1364
##    m0[12]  30.15  30.17 2.55 2.44   26.10  34.40 1.00     1954     1363
##    m0[13]  66.03  66.03 3.68 3.68   60.00  72.00 1.01     2038     1522
##    m0[14]  36.88  36.94 2.01 1.94   33.54  40.11 1.00     2301     1512
##    m0[15]   7.84   7.85 4.08 4.12    0.89  14.42 1.00     1388     1281
##    m0[16]  17.56  17.56 2.41 2.35   13.58  21.50 1.00     1835     1202
##    m0[17]  14.51  14.39 3.38 3.38    9.20  20.03 1.00     1938     1283
##    m0[18]  24.30  24.30 1.78 1.70   21.47  27.27 1.00     2044     1560
##    m0[19]  11.27  11.22 4.49 4.60    4.06  18.41 1.00     1929     1410
##    m0[20]  43.60  43.63 1.85 1.77   40.57  46.58 1.00     2144     1218
##    m0[21]  24.47  24.45 1.82 1.74   21.53  27.46 1.00     2038     1613
##    m0[22]  22.88  22.85 1.50 1.45   20.39  25.31 1.00     2104     1571
##    m0[23]  30.31  30.31 2.57 2.45   26.24  34.58 1.00     1953     1363
##    m0[24]  16.00  16.02 1.84 1.78   13.01  19.00 1.00     2011     1460
##    m0[25]  58.04  57.96 2.66 2.54   53.72  62.49 1.00     2125     1283
##    m0[26]  44.61  44.64 1.93 1.87   41.44  47.70 1.00     2149     1233
##    m0[27]  56.89  56.84 2.57 2.42   52.69  61.21 1.00     2120     1251
##    m0[28]  22.37  22.32 3.09 3.27   17.45  27.46 1.00     2030     1751
##    m0[29]  21.98  21.96 2.16 2.17   18.58  25.54 1.00     2891     1621
##    m0[30]  15.40  15.44 2.81 2.74   10.93  20.00 1.00     2890     1523
##    m0[31]  35.86  35.85 1.49 1.41   33.41  38.28 1.00     2296     1364
##    m0[32]  61.25  61.14 3.04 2.95   56.30  66.31 1.00     1969     1409
##    m0[33]  22.29  22.25 3.07 3.24   17.40  27.33 1.00     2032     1751
##    m0[34]   8.74   8.76 3.73 3.52    2.74  14.82 1.00     1612     1365
##    m0[35]  19.24  19.19 1.79 1.74   16.40  22.19 1.00     1941     1619
##    m0[36]  16.85  16.89 2.66 2.61   12.55  21.18 1.00     2890     1522
##    m0[37]  27.35  27.34 2.64 2.65   23.13  31.64 1.00     1872     1414
##    m0[38]  13.15  13.16 2.60 2.46    8.69  17.45 1.00     1721     1605
##    m0[39]  22.03  22.02 1.43 1.41   19.71  24.35 1.01     2101     1542
##    m0[40]  -8.90  -8.82 3.74 3.62  -15.08  -2.72 1.00     2322     1406
##    m0[41]  18.65  18.63 3.42 3.23   12.98  24.29 1.00     2315     1325
##    m0[42]  25.88  25.88 2.21 2.13   22.30  29.56 1.00     1925     1487
##    m0[43]  30.52  30.56 2.37 2.24   26.59  34.45 1.00     2282     1469
##    m0[44]  19.94  19.90 1.62 1.64   17.38  22.61 1.00     1963     1682
##    m0[45]  56.29  56.28 2.67 2.56   51.98  60.68 1.00     1988     1416
##    m0[46]  45.94  45.97 2.03 1.97   42.57  49.17 1.00     2166     1326
##    m0[47]  52.45  52.48 2.79 2.73   47.95  56.92 1.01     2051     1521
##    m0[48]  54.06  54.06 2.54 2.44   50.03  58.26 1.00     1998     1375
##    m0[49]  21.26  21.26 1.44 1.45   18.94  23.59 1.00     2046     1662
##    m0[50]   3.19   3.14 2.94 2.90   -1.52   8.07 1.00     1981     1464
##    y0[1]   32.05  32.05 5.71 5.60   22.58  41.35 1.00     1840     1920
##    y0[2]   39.87  39.96 4.74 4.72   31.90  47.42 1.00     2021     1618
##    y0[3]    5.92   5.90 5.24 5.15   -2.66  14.73 1.00     1793     1895
##    y0[4]    8.43   8.43 5.98 5.93   -1.30  18.13 1.00     1993     1710
##    y0[5]   57.86  57.82 5.47 5.26   49.06  67.06 1.00     1925     1744
##    y0[6]    2.83   2.74 5.65 5.61   -6.39  12.38 1.00     1974     1760
##    y0[7]   16.87  16.93 6.46 6.26    5.85  27.24 1.00     2042     1851
##    y0[8]   24.77  24.73 4.88 4.85   16.74  32.92 1.00     1517     1820
##    y0[9]   26.83  26.82 5.23 5.20   18.32  35.32 1.00     2029     1899
##    y0[10]  48.50  48.53 4.91 4.85   40.46  56.59 1.00     2025     2013
##    y0[11]   3.68   3.81 5.20 5.14   -5.10  12.10 1.00     1937     1729
##    y0[12]  30.17  30.10 5.28 5.25   21.39  38.97 1.00     1864     1720
##    y0[13]  65.93  65.83 5.76 5.65   56.41  75.70 1.00     1869     1639
##    y0[14]  36.86  36.85 5.08 5.18   28.63  45.18 1.00     2238     2099
##    y0[15]   7.81   7.98 5.94 5.76   -1.97  17.25 1.00     1725     1622
##    y0[16]  17.58  17.65 5.10 4.98    9.13  25.83 1.00     1850     1782
##    y0[17]  14.52  14.57 5.63 5.72    5.43  23.97 1.00     1866     1702
##    y0[18]  24.42  24.35 4.73 4.71   16.67  32.08 1.00     1488     1684
##    y0[19]  11.37  11.45 6.17 5.82    1.37  21.70 1.00     1973     1701
##    y0[20]  43.64  43.52 4.86 4.70   35.95  51.34 1.00     2087     1970
##    y0[21]  24.44  24.37 4.85 4.77   16.54  32.33 1.00     1843     1891
##    y0[22]  22.86  22.92 4.71 4.66   14.95  30.46 1.00     1917     1971
##    y0[23]  30.25  30.20 5.21 5.09   22.02  38.80 1.00     1965     1868
##    y0[24]  15.99  15.81 5.03 4.79    8.25  24.69 1.00     1650     1912
##    y0[25]  58.07  58.08 5.40 5.38   49.07  67.23 1.00     2117     1539
##    y0[26]  44.69  44.59 4.96 5.01   36.62  52.94 1.00     2033     1800
##    y0[27]  56.73  56.75 5.33 5.09   47.97  65.14 1.00     1994     1945
##    y0[28]  22.35  22.25 5.56 5.29   13.55  31.40 1.00     2009     1911
##    y0[29]  22.02  22.03 5.06 5.06   14.10  30.20 1.00     2136     1808
##    y0[30]  15.41  15.41 5.40 5.16    6.35  24.08 1.00     2181     1880
##    y0[31]  35.81  35.84 4.85 4.92   27.83  43.79 1.00     1883     1872
##    y0[32]  61.12  61.21 5.40 5.34   52.21  69.92 1.00     1838     1606
##    y0[33]  22.41  22.57 5.58 5.57   13.08  31.42 1.00     2058     1848
##    y0[34]   8.72   8.55 6.05 6.05   -0.99  18.98 1.00     1811     1729
##    y0[35]  19.06  19.08 4.90 4.82   10.98  27.15 1.00     2170     1895
##    y0[36]  16.90  16.82 5.26 5.11    7.84  25.43 1.00     2253     1879
##    y0[37]  27.20  27.18 5.22 5.04   18.80  36.18 1.00     2081     1921
##    y0[38]  13.15  13.17 5.23 5.17    4.54  21.73 1.00     2053     1754
##    y0[39]  22.03  22.10 4.79 4.64   14.35  29.89 1.00     1873     1949
##    y0[40]  -8.95  -9.00 5.86 5.72  -18.40   0.48 1.00     2182     1801
##    y0[41]  18.67  18.56 5.73 5.63    9.31  28.24 1.00     2071     1624
##    y0[42]  25.94  25.83 5.03 4.96   17.72  34.26 1.00     2074     1835
##    y0[43]  30.53  30.57 5.10 5.02   22.31  39.07 1.00     1991     1877
##    y0[44]  19.98  19.96 4.86 4.83   12.03  27.89 1.00     2031     1767
##    y0[45]  56.36  56.35 5.30 5.04   47.58  65.06 1.00     1862     1808
##    y0[46]  45.75  45.70 4.97 4.72   37.83  53.93 1.00     2176     2054
##    y0[47]  52.39  52.33 5.38 5.34   43.84  61.22 1.00     1978     1932
##    y0[48]  54.02  54.01 5.13 5.08   45.56  62.22 1.00     1993     1892
##    y0[49]  21.33  21.28 4.64 4.38   13.63  29.16 1.00     1920     1868
##    y0[50]   3.25   3.27 5.41 5.28   -5.61  12.52 1.00     2043     2099

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

estimate as class have relation

all b0l,b1l have a distribution  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)

ex12-2.stan

class have relation

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  real b00;
  real<lower=0,upper=100> sb0;
  vector[K] b0;
  real b10;
  real<lower=0,upper=100> sb1;
  vector[K] b1;
  real<lower=0,upper=100> s;
}
model{
  b0~normal(b00,sb0);
  b1~normal(b10,sb1);
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-2.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -314.534 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -36.5232      0.301266       99672.6       0.141        0.77      139    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199       121.716      0.391023   1.47123e+14     0.08681      0.5203      311    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      220       139.916      0.289332   4.96632e+15       1e-12       0.001      406  LS failed, Hessian reset  
## Finished in  0.1 seconds.
try(print(mle))
## Error : Fitting failed. Unable to print.
mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 3 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 3 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 4 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 4 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 4 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 1.3 seconds.
## Chain 2 finished in 1.3 seconds.
## Chain 4 finished in 1.3 seconds.
## Chain 3 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 3 finished in 1.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.5 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -123.72 -123.26 4.24 4.07 -131.31 -117.68 1.00     2474     4034
##    b00       8.45    8.46 3.40 3.04    2.99   13.78 1.00     7156     5013
##    sb0       8.62    8.02 3.35 2.76    4.55   14.78 1.00     4184     4309
##    b0[1]    10.44   10.43 3.96 3.88    3.99   17.11 1.00     8744     6284
##    b0[2]     5.03    5.03 4.31 4.22   -1.98   12.14 1.00     7276     5822
##    b0[3]    13.11   13.15 2.80 2.74    8.55   17.71 1.00     6169     5776
##    b0[4]    -5.80   -5.88 3.99 3.95  -12.22    0.86 1.00     6285     4906
##    b0[5]    13.44   13.35 4.13 4.10    6.84   20.35 1.00     7094     6068
##    b0[6]    14.48   14.47 3.51 3.47    8.73   20.25 1.00     6966     5527
##    b0[7]     7.40    7.42 3.61 3.47    1.47   13.38 1.00     8428     5679
##    b0[8]     6.74    6.83 4.32 4.17   -0.60   13.67 1.00     7687     6211
##    b0[9]    10.50   10.54 3.77 3.65    4.24   16.64 1.00     8060     6066
##    b10       1.67    1.68 0.43 0.38    0.98    2.36 1.00     6651     4829
##    sb1       1.18    1.10 0.43 0.32    0.68    1.97 1.00     6527     4576
##    b1[1]     2.74    2.75 0.29 0.29    2.26    3.22 1.00     8492     6055
##    b1[2]     0.95    0.94 0.32 0.31    0.43    1.47 1.00     7462     6178
##    b1[3]     2.17    2.17 0.25 0.24    1.77    2.58 1.00     7221     5744
##    b1[4]     1.80    1.81 0.27 0.26    1.35    2.23 1.00     6671     5405
##    b1[5]     2.71    2.72 0.30 0.30    2.20    3.21 1.00     7494     6070
##    b1[6]     2.20    2.20 0.24 0.23    1.81    2.59 1.00     7460     5214
##    b1[7]     1.56    1.56 0.34 0.34    1.00    2.12 1.00     8377     5988
##    b1[8]    -0.14   -0.15 0.37 0.36   -0.74    0.49 1.00     7696     5736
##    b1[9]     1.05    1.05 0.33 0.33    0.52    1.60 1.00     8216     6128
##    s         4.52    4.45 0.59 0.57    3.69    5.60 1.00     5091     5985
##    m0[1]    32.02   32.00 3.40 3.35   26.38   37.60 1.00     8126     6238
##    m0[2]    39.68   39.70 1.63 1.63   37.04   42.33 1.00     8875     5062
##    m0[3]     5.57    5.59 2.30 2.23    1.78    9.30 1.00     8821     6448
##    m0[4]     8.47    8.49 3.43 3.30    2.80   14.16 1.00     8413     5868
##    m0[5]    58.02   57.99 2.52 2.54   53.94   62.16 1.00     8305     6408
##    m0[6]     4.39    4.33 3.27 3.24   -0.93    9.87 1.00     8504     6460
##    m0[7]    16.23   16.14 3.87 3.82   10.05   22.72 1.00     7094     6361
##    m0[8]    25.25   25.27 1.84 1.80   22.24   28.28 1.00     8373     6156
##    m0[9]    27.60   27.59 2.37 2.32   23.78   31.50 1.00     8430     6514
##    m0[10]   48.44   48.44 2.24 2.23   44.79   52.14 1.00     8447     6094
##    m0[11]    5.88    5.83 2.60 2.54    1.71   10.17 1.00     6489     5230
##    m0[12]   29.56   29.59 2.57 2.51   25.30   33.79 1.00     8838     6296
##    m0[13]   65.12   65.18 3.72 3.65   58.97   71.22 1.00     8827     6172
##    m0[14]   36.40   36.43 1.93 1.91   33.23   39.56 1.00     7800     5755
##    m0[15]    6.29    6.37 3.37 3.25    0.58   11.72 1.00     7783     6676
##    m0[16]   17.58   17.57 2.30 2.27   13.81   21.36 1.00     8101     6628
##    m0[17]   13.21   13.23 3.00 2.90    8.28   18.10 1.00     8043     5919
##    m0[18]   24.62   24.65 1.73 1.69   21.80   27.47 1.00     8338     6055
##    m0[19]   11.75   11.72 3.86 3.78    5.45   18.23 1.00     8760     6389
##    m0[20]   43.50   43.51 1.86 1.83   40.48   46.54 1.00     8703     5623
##    m0[21]   24.81   24.84 1.76 1.72   21.93   27.70 1.00     8349     6228
##    m0[22]   22.97   22.96 1.51 1.46   20.49   25.40 1.00     8237     5028
##    m0[23]   29.70   29.73 2.59 2.53   25.42   33.94 1.00     8837     6310
##    m0[24]   16.94   16.90 1.91 1.86   13.82   20.10 1.00     7975     5092
##    m0[25]   58.21   58.18 2.54 2.55   54.10   62.37 1.00     8302     6173
##    m0[26]   44.52   44.53 1.93 1.92   41.39   47.67 1.00     8647     5684
##    m0[27]   57.02   57.01 2.45 2.45   53.04   61.02 1.00     8319     6214
##    m0[28]   22.71   22.65 3.08 3.00   17.78   27.74 1.00     8412     6197
##    m0[29]   21.57   21.58 2.04 2.02   18.23   24.96 1.00     6358     5809
##    m0[30]   14.90   14.94 2.62 2.58   10.63   19.22 1.00     6173     5631
##    m0[31]   35.65   35.67 1.51 1.50   33.20   38.10 1.00     8852     5134
##    m0[32]   61.02   61.01 2.86 2.83   56.36   65.69 1.00     8262     5445
##    m0[33]   22.63   22.57 3.06 2.99   17.73   27.64 1.00     8420     6097
##    m0[34]    9.02    9.04 3.23 3.17    3.78   14.40 1.00     7418     6338
##    m0[35]   18.72   18.74 1.70 1.65   15.91   21.46 1.00     7850     6647
##    m0[36]   16.37   16.40 2.49 2.43   12.35   20.45 1.00     6193     5630
##    m0[37]   27.31   27.28 2.70 2.67   22.88   31.74 1.00     8025     5571
##    m0[38]   13.45   13.46 2.36 2.31    9.62   17.35 1.00     7962     6548
##    m0[39]   21.97   21.97 1.45 1.39   19.57   24.28 1.00     8209     5373
##    m0[40]   -5.29   -5.36 3.92 3.88  -11.57    1.27 1.00     6289     4942
##    m0[41]   17.61   17.58 3.22 3.19   12.30   22.91 1.00     6942     5728
##    m0[42]   25.75   25.76 2.23 2.16   22.05   29.45 1.00     8854     5711
##    m0[43]   29.84   29.85 2.26 2.22   26.15   33.54 1.00     7298     5730
##    m0[44]   19.54   19.57 1.58 1.53   16.93   22.10 1.00     7910     6227
##    m0[45]   56.01   55.99 2.54 2.49   51.86   60.13 1.00     8286     5270
##    m0[46]   45.87   45.87 2.03 2.01   42.58   49.22 1.00     8577     5753
##    m0[47]   51.88   51.92 2.86 2.81   47.14   56.52 1.00     9043     5148
##    m0[48]   53.75   53.74 2.42 2.37   49.81   57.69 1.00     8283     5458
##    m0[49]   21.07   21.09 1.45 1.41   18.68   23.40 1.00     8210     5651
##    m0[50]    4.54    4.51 2.99 2.92   -0.37    9.52 1.00     8633     6172
##    y0[1]    31.97   31.91 5.66 5.55   22.79   41.23 1.00     7778     7766
##    y0[2]    39.77   39.76 4.87 4.76   31.73   47.89 1.00     8342     7861
##    y0[3]     5.61    5.57 5.07 4.92   -2.66   13.89 1.00     8037     7613
##    y0[4]     8.50    8.54 5.66 5.60   -0.79   17.62 1.00     7871     7393
##    y0[5]    58.03   58.02 5.18 5.10   49.55   66.53 1.00     8243     7442
##    y0[6]     4.39    4.31 5.66 5.57   -4.62   13.75 1.00     7883     7402
##    y0[7]    16.17   16.24 6.06 5.78    6.27   26.27 1.00     6878     7230
##    y0[8]    25.19   25.29 4.98 4.80   16.97   33.20 1.00     8209     7187
##    y0[9]    27.49   27.47 5.14 5.06   18.97   35.99 1.00     7292     7271
##    y0[10]   48.47   48.48 5.09 4.89   40.12   56.94 1.00     8290     7837
##    y0[11]    5.85    5.72 5.24 5.18   -2.63   14.52 1.00     7830     7396
##    y0[12]   29.52   29.45 5.25 5.09   20.91   38.23 1.00     8493     7789
##    y0[13]   65.17   65.26 5.90 5.73   55.52   74.73 1.00     8230     7849
##    y0[14]   36.35   36.40 4.99 4.85   28.14   44.55 1.00     7970     7757
##    y0[15]    6.34    6.37 5.64 5.54   -3.11   15.46 1.00     8018     7295
##    y0[16]   17.60   17.51 5.11 5.01    9.39   25.99 1.00     7777     7263
##    y0[17]   13.20   13.19 5.40 5.34    4.33   22.12 1.00     7789     7808
##    y0[18]   24.69   24.69 4.92 4.86   16.68   32.83 1.00     7981     8098
##    y0[19]   11.78   11.85 6.09 6.05    1.52   21.65 1.00     8346     7375
##    y0[20]   43.55   43.51 4.93 4.86   35.53   51.58 1.00     8138     7794
##    y0[21]   24.89   24.94 4.89 4.82   16.82   32.78 1.00     8370     7850
##    y0[22]   23.03   23.03 4.79 4.69   15.26   30.87 1.00     7878     7842
##    y0[23]   29.77   29.77 5.22 5.14   21.27   38.39 1.00     8240     8098
##    y0[24]   16.89   16.83 4.96 4.91    8.78   24.99 1.00     8019     7504
##    y0[25]   58.24   58.27 5.23 5.12   49.66   66.84 1.00     8108     7101
##    y0[26]   44.67   44.68 4.82 4.71   36.84   52.53 1.00     7954     7626
##    y0[27]   57.02   57.05 5.17 5.12   48.60   65.50 1.00     7939     7647
##    y0[28]   22.72   22.72 5.50 5.52   13.70   31.73 1.00     7759     7674
##    y0[29]   21.59   21.69 5.03 4.92   13.31   29.63 1.00     7534     7540
##    y0[30]   14.93   15.04 5.25 5.17    6.12   23.40 1.00     7651     7808
##    y0[31]   35.68   35.66 4.79 4.75   27.82   43.48 1.00     8558     7674
##    y0[32]   61.03   60.99 5.34 5.21   52.29   69.81 1.00     7870     7778
##    y0[33]   22.67   22.65 5.47 5.41   13.83   31.72 1.00     8005     7076
##    y0[34]    8.96    8.86 5.55 5.50   -0.10   18.10 1.00     7593     7594
##    y0[35]   18.69   18.73 4.92 4.74   10.53   26.70 1.00     7900     7498
##    y0[36]   16.35   16.36 5.23 5.16    7.65   24.92 1.00     7336     7398
##    y0[37]   27.31   27.26 5.32 5.23   18.60   36.15 1.00     7784     7633
##    y0[38]   13.43   13.39 5.14 5.04    5.08   22.02 1.00     8160     7544
##    y0[39]   21.97   21.97 4.77 4.62   14.13   29.86 1.00     7064     7320
##    y0[40]   -5.23   -5.26 6.06 5.91  -15.09    4.83 1.00     7198     6895
##    y0[41]   17.55   17.57 5.61 5.48    8.39   26.85 1.00     8214     7125
##    y0[42]   25.72   25.67 5.10 5.03   17.53   34.16 1.00     8198     7550
##    y0[43]   29.85   29.83 5.10 5.06   21.52   38.16 1.00     7791     7972
##    y0[44]   19.56   19.60 4.80 4.67   11.62   27.42 1.00     8096     7875
##    y0[45]   55.97   56.09 5.21 5.11   47.45   64.32 1.00     8078     7710
##    y0[46]   45.84   45.77 4.94 4.88   37.78   54.10 1.00     7572     7481
##    y0[47]   51.85   51.88 5.41 5.22   43.04   60.84 1.00     8736     7203
##    y0[48]   53.68   53.74 5.13 5.03   45.32   62.07 1.00     8136     6479
##    y0[49]   21.05   21.00 4.76 4.66   13.27   28.92 1.00     7153     7071
##    y0[50]    4.62    4.49 5.47 5.37   -4.18   13.72 1.00     8256     7782

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

generalized linear mixed model

(X,y)i=1-n
b[b0,b1,...]
linear model    y~N(Xb,s)  
generalized linear model    y~dist.(m=link(Xb),s)  

fixed effect    b0, b1,...  
individual random effect   b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...  
class c=1-k  
class effect    b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...  

for y=dist.(m,s)
random intercept model    m=b0+r0i+b1*x, m=b0+r0c+b1*x  
random coefficient model  m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x  
mixed model   m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x  

note  
@ yi=b0+b1*xi+r0i is not useful, ri is included in s  
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
  but variance is larger than mean (over dispersion),
  random intercept model is useful
n=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)

ex13-0.stan

generalized linear model, poisson regression

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~poisson_log(b0+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-0.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -8.45564e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       21       1184.69   0.000250358     0.0562033      0.9629      0.9629       61    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    1184.69
##    b0         2.80
##    b1         0.08
##    m0[1]      3.19
##    m0[2]      3.26
##    m0[3]      3.56
##    m0[4]      3.32
##    m0[5]      3.31
##    m0[6]      3.43
##    m0[7]      3.30
##    m0[8]      3.11
##    m0[9]      3.40
##    m0[10]     2.82
##    m0[11]     2.81
##    m0[12]     3.46
##    m0[13]     3.46
##    m0[14]     2.98
##    m0[15]     2.86
##    m0[16]     3.28
##    m0[17]     3.15
##    m0[18]     3.52
##    m0[19]     3.50
##    m0[20]     2.86
##    y0[1]     19.00
##    y0[2]     23.00
##    y0[3]     31.00
##    y0[4]     24.00
##    y0[5]     34.00
##    y0[6]     30.00
##    y0[7]     32.00
##    y0[8]     20.00
##    y0[9]     38.00
##    y0[10]    13.00
##    y0[11]    17.00
##    y0[12]    29.00
##    y0[13]    33.00
##    y0[14]    21.00
##    y0[15]    24.00
##    y0[16]    21.00
##    y0[17]    32.00
##    y0[18]    34.00
##    y0[19]    34.00
##    y0[20]    17.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   1183.21 1183.95 4.12 0.77 1181.28 1184.64 1.01      243      148
##    b0        2.78    2.79 0.15 0.11    2.60    2.96 1.01      276      122
##    b1        0.08    0.08 0.02 0.02    0.05    0.11 1.01      297      138
##    m0[1]     3.18    3.19 0.06 0.05    3.10    3.26 1.01      374      147
##    m0[2]     3.25    3.25 0.05 0.05    3.17    3.33 1.01      461      159
##    m0[3]     3.56    3.56 0.08 0.07    3.43    3.68 1.01      495      266
##    m0[4]     3.32    3.32 0.05 0.05    3.24    3.39 1.01      671      243
##    m0[5]     3.30    3.30 0.05 0.05    3.22    3.38 1.01      594      238
##    m0[6]     3.42    3.43 0.05 0.05    3.33    3.51 1.00     1643     1532
##    m0[7]     3.30    3.30 0.05 0.05    3.22    3.37 1.01      570      236
##    m0[8]     3.11    3.11 0.08 0.06    3.01    3.20 1.01      319      127
##    m0[9]     3.40    3.40 0.05 0.05    3.31    3.48 1.00     1873     1457
##    m0[10]    2.80    2.81 0.15 0.11    2.62    2.97 1.01      276      120
##    m0[11]    2.79    2.81 0.15 0.11    2.61    2.97 1.01      276      120
##    m0[12]    3.46    3.46 0.06 0.06    3.36    3.56 1.00      982      666
##    m0[13]    3.46    3.46 0.06 0.06    3.36    3.56 1.00      974      666
##    m0[14]    2.97    2.98 0.11 0.08    2.85    3.10 1.01      286      123
##    m0[15]    2.84    2.85 0.14 0.10    2.67    3.00 1.01      277      123
##    m0[16]    3.27    3.27 0.05 0.05    3.19    3.35 1.01      503      179
##    m0[17]    3.14    3.15 0.07 0.05    3.05    3.23 1.01      340      116
##    m0[18]    3.52    3.52 0.07 0.06    3.40    3.62 1.01      580      322
##    m0[19]    3.50    3.50 0.07 0.06    3.38    3.60 1.01      655      391
##    m0[20]    2.85    2.86 0.14 0.10    2.68    3.01 1.01      277      121
##    y0[1]    24.02   24.00 5.21 5.93   16.00   33.00 1.00     1493     1602
##    y0[2]    25.83   26.00 5.24 5.93   17.00   35.00 1.00     1743     1864
##    y0[3]    35.38   35.00 6.50 5.93   25.00   47.00 1.00     1541     1813
##    y0[4]    27.68   27.00 5.30 5.93   19.00   36.05 1.00     1906     2038
##    y0[5]    27.13   27.00 5.60 5.93   18.00   36.00 1.00     1691     1828
##    y0[6]    30.72   30.50 5.71 5.19   22.00   41.00 1.00     1968     1722
##    y0[7]    27.18   27.00 5.36 5.93   19.00   36.00 1.00     1876     1885
##    y0[8]    22.36   22.00 5.05 4.45   14.00   31.00 1.01     1015     1039
##    y0[9]    29.89   30.00 5.64 5.93   21.00   40.00 1.00     1969     1858
##    y0[10]   16.59   16.00 4.61 4.45    9.00   24.00 1.00      650      463
##    y0[11]   16.39   16.00 4.54 4.45    9.00   24.00 1.00      630      358
##    y0[12]   31.80   32.00 6.02 5.93   22.00   42.00 1.00     2125     1941
##    y0[13]   31.87   32.00 5.79 5.93   23.00   42.00 1.00     1996     1941
##    y0[14]   19.52   19.00 4.83 4.45   12.00   28.00 1.00      846      511
##    y0[15]   17.25   17.00 4.71 4.45   10.00   25.00 1.00      654      493
##    y0[16]   26.39   26.00 5.38 5.93   18.00   36.00 1.00     1722     1847
##    y0[17]   23.11   23.00 5.04 4.45   15.00   31.00 1.00     1185      923
##    y0[18]   33.70   34.00 6.26 5.93   24.00   44.00 1.00     1955     1254
##    y0[19]   32.93   33.00 6.06 5.93   23.00   43.00 1.00     1276     1800
##    y0[20]   17.19   17.00 4.52 4.45   10.00   25.00 1.00      650      401

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

ex13-1.stan

generalized linear mixed model, poisson regression

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
  real<lower=0> sr0;
  vector[N] r0;
}
model{
  r0~normal(0,sr0);
  for(i in 1:N)
    y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+r0[i]+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-1.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -7.77186e+07 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       99       23721.8      0.404198       109.714      0.9746      0.9746      153    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199       24031.5      0.425447   1.22229e+08      0.1516           1      314    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      299       24418.7       0.80365    9.7267e+18       1.741      0.5183      481    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      399       24678.1     0.0183169   1.03074e+29     0.09589     0.09589      669    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      499       24817.7        0.2255   3.83185e+39       2.326      0.3774      857    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      541       24829.3      0.044588   2.15164e+41      0.8437      0.8437      923    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__   24829.30
##    b0         2.51
##    b1         0.17
##    sr0        0.00
##    r0[1]      0.00
##    r0[2]      0.00
##    r0[3]      0.00
##    r0[4]      0.00
##    r0[5]      0.00
##    r0[6]      0.00
##    r0[7]      0.00
##    r0[8]      0.00
##    r0[9]      0.00
##    r0[10]     0.00
##    r0[11]     0.00
##    r0[12]     0.00
##    r0[13]     0.00
##    r0[14]     0.00
##    r0[15]     0.00
##    r0[16]     0.00
##    r0[17]     0.00
##    r0[18]     0.00
##    r0[19]     0.00
##    r0[20]     0.00
##    m0[1]      3.36
##    m0[2]      3.51
##    m0[3]      4.17
##    m0[4]      3.65
##    m0[5]      3.62
##    m0[6]      3.87
##    m0[7]      3.61
##    m0[8]      3.20
##    m0[9]      3.81
##    m0[10]     2.56
##    m0[11]     2.54
##    m0[12]     3.95
##    m0[13]     3.95
##    m0[14]     2.92
##    m0[15]     2.64
##    m0[16]     3.55
##    m0[17]     3.28
##    m0[18]     4.07
##    m0[19]     4.02
##    m0[20]     2.65
##    y0[1]     24.00
##    y0[2]     32.00
##    y0[3]     63.00
##    y0[4]     34.00
##    y0[5]     41.00
##    y0[6]     41.00
##    y0[7]     40.00
##    y0[8]     23.00
##    y0[9]     58.00
##    y0[10]    14.00
##    y0[11]    13.00
##    y0[12]    50.00
##    y0[13]    52.00
##    y0[14]    17.00
##    y0[15]     8.00
##    y0[16]    30.00
##    y0[17]    30.00
##    y0[18]    68.00
##    y0[19]    61.00
##    y0[20]     9.00
mcmc=goMCMC(mdl,data,wrm=500)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 1.2 seconds.
## Chain 2 finished in 1.1 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 finished in 1.3 seconds.
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 finished in 1.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.2 seconds.
## Total execution time: 1.5 seconds.
seeMCMC(mcmc,'[m,r]')
##  variable     mean   median    sd   mad       q5      q95 rhat ess_bulk
##    lp__   23773.49 23771.70 14.54 14.23 23752.20 23799.80 1.08       51
##    b0         2.80     2.79  0.02  0.02     2.75     2.83 1.01      768
##    b1         0.08     0.08  0.00  0.00     0.07     0.08 1.01      836
##    sr0        0.01     0.01  0.01  0.01     0.00     0.02 1.07       59
##    r0[1]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3212
##    r0[2]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2708
##    r0[3]      0.00     0.00  0.01  0.01    -0.02     0.02 1.04     2613
##    r0[4]      0.00     0.00  0.01  0.01    -0.02     0.02 1.01     3172
##    r0[5]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2719
##    r0[6]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2865
##    r0[7]      0.00     0.00  0.01  0.01    -0.02     0.02 1.03     2699
##    r0[8]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2589
##    r0[9]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2598
##    r0[10]     0.00     0.00  0.01  0.01    -0.02     0.02 1.03     2887
##    r0[11]     0.00     0.00  0.01  0.01    -0.02     0.02 1.03     3730
##    r0[12]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2705
##    r0[13]     0.00     0.00  0.01  0.01    -0.02     0.02 1.04     2441
##    r0[14]     0.00     0.00  0.01  0.01    -0.02     0.02 1.04     3104
##    r0[15]     0.00     0.00  0.01  0.01    -0.02     0.02 1.03     3611
##    r0[16]     0.00     0.00  0.01  0.01    -0.02     0.02 1.03     3203
##    r0[17]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3015
##    r0[18]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3160
##    r0[19]     0.00     0.00  0.01  0.01    -0.02     0.02 1.03     2525
##    r0[20]     0.00     0.00  0.01  0.01    -0.02     0.02 1.04     2418
##    m0[1]      3.19     3.19  0.02  0.01     3.17     3.21 1.00     1524
##    m0[2]      3.26     3.25  0.02  0.01     3.23     3.28 1.00     2167
##    m0[3]      3.56     3.56  0.02  0.02     3.53     3.59 1.00     1573
##    m0[4]      3.32     3.32  0.01  0.01     3.30     3.34 1.00     2428
##    m0[5]      3.31     3.31  0.01  0.01     3.28     3.33 1.00     2210
##    m0[6]      3.43     3.43  0.02  0.02     3.40     3.45 1.00     2419
##    m0[7]      3.30     3.30  0.02  0.01     3.28     3.33 1.01     2458
##    m0[8]      3.11     3.11  0.02  0.02     3.09     3.14 1.01     1186
##    m0[9]      3.40     3.40  0.02  0.01     3.37     3.42 1.00     2258
##    m0[10]     2.82     2.82  0.03  0.02     2.77     2.86 1.01      858
##    m0[11]     2.81     2.81  0.03  0.03     2.76     2.85 1.01      929
##    m0[12]     3.46     3.46  0.02  0.02     3.43     3.49 1.00     1934
##    m0[13]     3.46     3.46  0.02  0.02     3.43     3.49 1.00     1939
##    m0[14]     2.98     2.98  0.02  0.02     2.95     3.01 1.01      995
##    m0[15]     2.86     2.85  0.02  0.02     2.81     2.89 1.01      928
##    m0[16]     3.28     3.28  0.02  0.01     3.25     3.30 1.00     2201
##    m0[17]     3.15     3.15  0.02  0.02     3.12     3.17 1.01     1351
##    m0[18]     3.51     3.51  0.02  0.02     3.49     3.54 1.00     2035
##    m0[19]     3.50     3.50  0.02  0.02     3.47     3.53 1.00     2359
##    m0[20]     2.86     2.86  0.02  0.02     2.82     2.90 1.00      901
##    y0[1]     24.13    24.00  4.94  4.45    16.00    32.00 1.00     1976
##    y0[2]     25.90    26.00  5.04  4.45    18.00    34.00 1.00     2003
##    y0[3]     35.09    35.00  5.83  5.93    26.00    45.00 1.00     1867
##    y0[4]     27.59    28.00  5.12  5.19    20.00    36.00 1.00     2074
##    y0[5]     27.43    27.00  5.50  5.93    19.00    37.00 1.00     2060
##    y0[6]     30.59    30.00  5.46  5.93    22.00    40.00 1.00     2094
##    y0[7]     26.99    27.00  5.15  4.45    19.00    36.00 1.00     1969
##    y0[8]     22.53    22.00  4.71  4.45    15.00    31.00 1.00     1866
##    y0[9]     29.72    29.00  5.56  5.93    21.00    39.00 1.00     2026
##    y0[10]    16.84    17.00  4.16  4.45    10.00    24.00 1.00     2103
##    y0[11]    16.65    16.00  4.08  4.45    10.00    24.00 1.00     1837
##    y0[12]    31.83    32.00  5.68  5.93    23.00    42.00 1.00     1998
##    y0[13]    32.10    32.00  5.76  5.93    23.00    42.00 1.00     2100
##    y0[14]    19.59    20.00  4.45  4.45    12.00    27.00 1.00     1967
##    y0[15]    17.32    17.00  4.27  4.45    11.00    25.00 1.00     2009
##    y0[16]    26.48    26.00  5.05  4.45    18.00    35.00 1.00     1944
##    y0[17]    23.39    23.00  4.94  4.45    15.00    32.00 1.00     2131
##    y0[18]    33.44    33.00  5.83  5.93    24.00    43.00 1.00     1985
##    y0[19]    32.86    33.00  5.73  5.93    24.00    43.00 1.00     1867
##    y0[20]    17.49    17.00  4.24  4.45    11.00    25.00 1.00     2022
##  ess_tail
##        88
##       284
##      1029
##       102
##      1053
##       728
##      1005
##       797
##       867
##      1070
##       962
##       989
##       552
##       693
##       944
##       915
##       618
##       808
##       919
##       876
##       656
##      1137
##       702
##       845
##      1259
##      1313
##      1539
##      1202
##      1523
##      1358
##      1303
##      1442
##      1053
##       280
##       303
##      1358
##      1222
##       605
##       289
##      1373
##      1540
##      1467
##      1408
##       291
##      1955
##      1995
##      1768
##      1979
##      2032
##      1882
##      1877
##      1855
##      1954
##      1926
##      1829
##      2019
##      2055
##      2067
##      1908
##      1904
##      1801
##      1961
##      1955
##      1954

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

categorical dist. with parameters following Dirichlet dist.

h~Dir(a), h[0,1], sum(h)=1, a[0,1], sum(a)=1 
y~Cat(h), y[0,1], sum(y)=1

ex15-1.stan

data {
  int N;
  int K;
  array[N] int<lower=1,upper=K> y;
}
parameters {
  simplex[K] a;
  simplex[K] h;
}
model {
  h~dirichlet(a);
  y~categorical(h);
}
library(gtools)
n=20
h=rdirichlet(n,c(0.5,0.3,0.2))
summary(h)
##        V1              V2              V3       
##  Min.   :0.008   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.130   1st Qu.:0.010   1st Qu.:0.004  
##  Median :0.413   Median :0.144   Median :0.020  
##  Mean   :0.449   Mean   :0.345   Mean   :0.206  
##  3rd Qu.:0.777   3rd Qu.:0.751   3rd Qu.:0.318  
##  Max.   :0.996   Max.   :0.985   Max.   :0.992
c0=c(1,2,3)
for(i in 1:n) y[i]=sample(c0,1,prob=h[i,])
table(y) |> prop.table()
## y
##    1    2    3 
## 0.50 0.35 0.15
data=list(N=n,K=ncol(h),y=y)

mdl=cmdstan_model('./ex15-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -24.9423 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        9      -20.4197   7.57839e-05   0.000195176      0.6424      0.6424       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -20.42
##      a[1]     0.39
##      a[2]     0.35
##      a[3]     0.26
##      h[1]     0.52
##      h[2]     0.35
##      h[3]     0.13
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -29.50 -29.16 1.45 1.23 -32.29 -27.75 1.00      953     1423
##      a[1]   0.36   0.34 0.18 0.20   0.10   0.68 1.00     1777     1354
##      a[2]   0.34   0.32 0.18 0.19   0.08   0.65 1.01     1369     1235
##      a[3]   0.30   0.28 0.16 0.17   0.07   0.60 1.00     1606     1421
##      h[1]   0.49   0.49 0.11 0.11   0.31   0.66 1.00     3231     1478
##      h[2]   0.35   0.34 0.10 0.10   0.19   0.52 1.00     2680     1702
##      h[3]   0.16   0.15 0.08 0.08   0.05   0.31 1.00     2122     1328

beta binomial dist.

y~betaB(n,a,b): y~B(n,p), p~beta(a,b)

ex15-2.stan

data {
  int N;
  int n;
  array[N] int y;
}
parameters {
  real<lower=0> a;
  real<lower=0> b;
}
model {
  a~cauchy(0,5);
  b~cauchy(0,5);
  y~beta_binomial(n,a,b);
}
generated quantities{
  int y1;
  y1=beta_binomial_rng(n,a,b);
}
n=30
a=3
b=4
p=rbeta(n,a,b)
n0=10
y=rbinom(n,n0,p)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.25    4.00    3.87    5.00    8.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex15-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -204.579 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -199.196    0.00183778    0.00380116      0.9979      0.9979        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -199.20
##      a        3.25
##      b        5.10
##      y1       1.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -197.07 -196.79 1.02 0.76 -198.91 -196.10 1.00      449      782
##      a       5.60    4.86 3.28 2.21    2.37   11.03 1.01      493      543
##      b       8.88    7.71 5.19 3.50    3.60   17.53 1.01      504      550
##      y1      3.89    4.00 2.00 1.48    1.00    7.00 1.00     1755     1684

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')